Useful R code 2

 

Date: 2020-11-13
R version: 3.5.0
*Corresponding author: matthew.malishev@gmail.com
This document can be found at https://github.com/darwinanddavis/UsefulCode

Overview

Same deal as Useful Code, but the second instalment because the first one has too much stuff in it and now runs slow.

Colour palettes

Colorspace

require(colorspace)
hcl_palettes(plot = TRUE)  # show all palettes

# https://cran.r-project.org/web/packages/colorspace/vignettes/colorspace.html
require(colorspace)
q4 <- qualitative_hcl(4, palette = "Dark 3")  # discrete
s9 <- sequential_hcl(9, "Purples 3")  # continuous
# for ggplot
scale_color_discrete_sequential(palette = "Purples 3", nmax = 6, order = 2:6)
# for colospace functions: hcl_palettes() %>% str hcl_palettes()['type']

Neon colour palettes

# https://www.shutterstock.com/blog/neon-color-palettes
neon1 <- c("#3B27BA", "#FF61BE", "#13CA91", "#FF9472")
neon2 <- c("#FFDEF3", "#FF61BE", "#3B55CE", "#35212A")
neon3 <- c("#FEA0FE", "#F85125", "#02B8A2", "#535EEB")
neon4 <- c("#535EEB", "#001437", "#C6BDEA", "#FFAA01")
scales::show_col(c(neon1, neon2, neon3, neon4))

Hexadecimal color code for transparency
See https://gist.github.com/lopspower/03fb1cc0ac9f32ef38f4.

require(colorspace)
require(stringr)
colv <- c("#004616", sequential_hcl(5, "Lajolla"))
str_sub(colv, 0, 1) <- "#66"  # add alpha opac to col vector

D3

Links
-
-
-

D3 and leaflet

# devtools::install_github('jcheng5/d3scatter')
require(pacman)
p_load(d3scatter, crosstalk, leaflet, tibble, httpuv)

sd <- SharedData$new(quakes[sample(nrow(quakes), 100), ])
# sd$data() %>% head

bscols(widths = c(12, 6, 6), filter_slider("stations", "Stations", sd, ~stations), leaflet(sd, width = "100%", 
    height = 400) %>% addTiles() %>% addCircleMarkers(lng = sd$data()[, "long"], lat = sd$data()[, "lat"], 
    stroke = F, fill = T, color = "red", fillOpacity = 0.5, radius = ~mag + 2, label = ~paste0("Depth: ", 
        as.character(depth))), d3scatter(sd, width = "100%", height = 400, ~mag, ~depth, color = ~stations))

Stations

Convert R code to D3 https://rstudio.github.io/r2d3/articles/visualization_options.html

Create calendar plot

# https://rstudio.github.io/r2d3/articles/gallery/calendar/ install.packages('r2d3')
require(r2d3)
require(readr)
require(dplyr)
require(colorspace)
require(scales)
require(stringr)

# col pal
col <- "PuBuGn"  # seq
col2 <- "Tropic"  # diverge

# seq
pal <- sequential_hcl(12, col)
# pal %>% show_col(borders = NA,labels=F)
paste0("\"", pal, "\"") %>% cat(sep = ",")

"#004533","#005C4E","#00726F","#008795","#0095B5","#56A0C8","#8FACD6","#B6BAE0","#D2CAE7","#E6DAEE","#F5EAF5","#FFF7FD"

# diverge
pal <- diverge_hcl(12, col2)
# pal %>% show_col(borders = NA,labels=F) paste0(''',pal,''') %>% cat(sep=',')

cal <- read_csv("https://raw.githubusercontent.com/rstudio/r2d3/master/vignettes/gallery/calendar/dji-latest.csv")

r2d3(data = cal, d3_version = 4, container = "div", options = list(start = 2006, end = 2011), script = "calendar.js")

Data frames

Reversing order of rows in dataframe

# df = data.frame
require(tidyverse)
df %>% map_df(rev)

Visualise data structure as tree

# explore package
require(DataExplorer)
require(palmerpenguins)
p <- penguins
plot_str(p)

dplyr basics

require(dplyr, gapminder)
pacman::p_load(gapminder)

# mutate
africa_ranked <- mutate(gapminder, African = continent == "Africa", RankPop = rank(desc(pop)))

africa_ranked %>% glimpse()  # visualise the data
africa_ranked %>% filter(continent == "Africa") %>% glimpse  # visualise just africa

# summarise data into one line
gapminder %>% summarise(MinYear = min(year, na.rm = T), MaxYear = max(year), CountryCount = n_distinct(country), 
    Counts = n())

gapminder %>% summarise(median(lifeExp))

# filter
require(gapminder)
gapminder %>% filter(continent == "Africa")

# group by
gapminder %>% group_by(continent) %>% summarise(median(lifeExp))


# group by continent and filter by year
life_cont_1992 <- gapminder %>% group_by(Continent = continent) %>% filter(year == 1992) %>% summarise(LifeExpect = median(lifeExp))

Execute unfriendly pipe functions inline in pipes

require(palmerpenguins)
require(dplyr)
p <- penguins

# %T>%
p %T>% glimpse %>% select(island)

# with()
p %>% with(lm(body_mass_g ~ flipper_length_mm)) %>% summary()

# %$% when var on lhs is undefined
require(magrittr)
data.frame(z = rnorm(100)) %$% ts.plot(z)

Generic functions

Convert character class to numeric (ideal when creating colour palettes to turn string cols in df to numeric)

require(dplyr)
set.seed(12)
df <- data.frame(X = LETTERS[sample(20)])
int_vec <- df$X %>% unlist %>% as.factor %>% as.integer  # converts to numbers 
int_vec
df$I <- int_vec
df

Pipe vector to multiple arguments

require(dplyr)

# as list
Sys.time() %>% list(format(., c("%y-%m", "%Y-%m", "%Y-%m")))

# use curly braces to keep original class
Sys.time() %>% {
    format(., c("%y-%m", "%Y-%m", "%Y-%m"))
}

Merge two data frames based on value and retain original number of rows

merge(a, b, by = "ID", sort = F)

HTML code

code <- "<!DOCTYPE html>
  <html>
<body>

<h1>My First Heading</h1>

<p>My first paragraph.</p>

</body>
</html>"

code <- paste(as.character(code), collapse = "\n")

write.table(code, file = "/Users/code.html", quote = FALSE, col.names = FALSE, row.names = FALSE)

Interactive plots

Clickme, NVD3, Polychart, rCharts, Rickshaw, and xCharts in R.
Link to collated Github page.

Leaflet

Interactive label options and custom tiles

require(leaflet)
require(dplyr)
require(geosphere)
require(htmltools)

setview <- c(7.369722, 12.354722)
mp <- data.frame(name = c("Melbourne", "Atlanta"), lat = c(-37.813629, 33.748997), lon = c(144.963058, 
    -84.387985))
latlon_matrix <- matrix(c(mp[, "lon"], mp[, "lat"]), ncol = 2)
custom_tile <- "http://a.sm.mapstack.stamen.com/(positron,(mapbox-water,$776699[hsl-color]),(buildings,$002bff[hsl-color]),(parks,$6abb9d[hsl-color]))/{z}/{x}/{y}.png"
colv <- "#4C3661"
opac <- 0.5
site_names <- mp$name
ttl <- "Debunking Flat Earth theory 101"
weblink <- "https://github.com/darwinanddavis"  # weblink
webname <- "My github"
href <- paste0("<b><a href=", weblink, ">", webname, "</a></b>")
text_label <- paste(sep = "<br/>", href, "606 5th Ave. S", "Seattle, WA 98138")
# label options
marker_label_opt <- labelOptions(textsize = "20px", opacity = 0.5, offset = c(0, 0))
text_label_opt <- labelOptions(noHide = T, direction = "top", textOnly = T, opacity = 1, offset = c(0, 
    0))

# title
tag.map.title <- tags$style(HTML(".leaflet-control.map-title { 
       transform: translate(-50%,20%);
       position: fixed !important;
       left: 50%;
       text-align: center;
       padding-left: 10px; 
       padding-right: 10px; 
       background: white; opacity: 0.7;
       font-weight: bold;
       font-size: 25px;
       }"))
title <- tags$div(tag.map.title, HTML(ttl))


# map
map <- gcIntermediate(latlon_matrix[1, ], latlon_matrix[2, ], n = 100, addStartEnd = T, sp = T) %>% leaflet() %>% 
    setView(setview[2], setview[1], zoom = 3) %>% addTiles(custom_tile) %>% addCircleMarkers(mp[, "lon"], 
    mp[, "lat"], radius = 10, stroke = TRUE, weight = 3, opacity = opac, color = colv, fillColor = colv, 
    label = paste(site_names), labelOptions = marker_label_opt) %>% addPolylines(color = colv, opacity = opac) %>% 
    addPopups(-122.327298, 47.597131, text_label, options = popupOptions(closeButton = FALSE, textOnly = T)) %>% 
    addLabelOnlyMarkers(setview[2], setview[1], label = text_label, labelOptions = text_label_opt) %>% 
    addControl("@darwinanddavis", position = "topright") %>% addControl(title, position = "topleft", 
    className = "map-title")
map

Lists

Transpose list (flip list elements)

l <- list(1:2, 3:4, 5:7, 8:10)
l
b <- data.table::transpose(l)
b

lengths for getting length of list indices

require(dplyr)
ls = list(rep(list(sample(50, replace = T)), 5))
ls %>% length
ls %>% lengths
lapply(ls, lengths)

Split list into smaller sublists

la = rep(list(1:5), 6)
names(la) = rep(LETTERS[1:3], 2)
u <- length(unique(names(la)))
n <- length(la)/u
split(la, rep(1:n, each = u))

# for when list has two elements in the name that change create a list of 10 letters with 5 lists in
# each
big_list <- rep(list(1:10), 5) %>% pmap(list)
names(big_list) <- LETTERS[1:10]
# to index the upper list
big_list["B"]  # 1
pluck(big_list, "B")  # 2
# to index the sublists
map(big_list["B"], 3)  # 1
bb_final <- list()  # 2
for (i in 1:10) {
    bb <- big_list["B"]
    bb_final <- c(bb_final, bb)
}
bb_final

Fill list elements with NAs to match length of longest element

# https://stackoverflow.com/questions/34570860/add-nas-to-make-all-list-elements-equal-length

# for single index list
set.seed(1)
ls = replicate(5, sample(1:100, 10), simplify = FALSE)
names(ls) = LETTERS[1:length(ls)]
lapply(ls, `length<-`, max(lengths(ls)))

# for sublists
ls = list(replicate(5, sample(1:100, 10), simplify = FALSE))
n.ticks = 20
fillvec = function(x) {
    nv = lapply(x, `length<-`, n.ticks)  # fill remaining vec with NAs to match total length
    rapply(nv, f = function(x) ifelse(is.na(x), 0, x), how = "replace")  # replace NAs with 0s
}
lapply(ls, fillvec)  # apply fillvec to list

Access list elements in loop by name/string

set.seed(12)

# inputs
time <- 5
time_vec <- 1:10
a_vec <- runif(10)
beta1_vec <- 1:10
beta2_vec <- 11:20
param_vec <- list(a_vec,beta1_vec,beta2_vec)
names(param_vec) <-c("alpha","beta1","beta2") 
params <- sapply(rep(NA,length(param_vec)),list) # create empty final params vector
names(params) <- names(param_vec)

# select parameter to test 
param_input <- "alpha" #beta1 #beta2

# run from here -----------------------------------------------------------
for(time in time_vec) {
  p_in = param_vec[`param_input`][[1]][time] # get parameter value by name
  # create new list of with updated param_input value
  params <- c(param_vec[-which(names(param_vec)==param_input)], # everything but param_input
              param_input = p_in # param_input
              )
  # get just the latest value
  # remove this if you want all list elements
  params <- sapply(params,function(x) x[1]) %>% as.numeric 
  # rename this new list
  names(params) <- c(names(param_vec)[-which(names(param_vec)==param_input)], # everything but param_input
                     param_input
                     )
  print(params)
} # end loop
params # each list element changes depending on user input 

Apply function to nested lists

ls = list(replicate(5, sample(1:100, 10), simplify = FALSE))
ls %>% glimpse
lapply(ls, lapply, mean)
lapply(ls, sapply, mean)  # return as one list 
rapply(ls, mean, how = "unlist")  # unlist, replace, or list

Loading packages

pacman

require(pacman)
p_load(dplyr, mapdeck)

Plotting

Hand drawn plotting using roughViz.js. Link to package page.

# install.packages('remotes') remotes::install_github('XiangyunHuang/roughviz')
require(roughviz)

Read in data

Read in csv data sources directly from web

# link to raw csv link on e.g. github
require(readr)
url <- "https://raw.githubusercontent.com/plotly/datasets/master/2011_february_aa_flight_paths.csv"
flights <- read_csv(url)

Regex

resource_type <- "algae"
# this regex expression
list.files(pattern = paste0("^", resource_type, "_[0,5]{1}_[0-9]{1,2}_hostpop50_predpop", "[0-9]{1,3}_rep[1-5]{1}\\.R$"))
# returns this begins with resource_type, either 0 or 5 as one integer, 0 to 9 as either one or two
# integers, 0 to 9 as one to three integers, and 1 to 5 as one integer
"algae_0_5_hostpop50_predpop5_rep1.R"
"algae_5_20_hostpop50_predpop30_rep2.R"
"algae_0_15_hostpop50_predpop150_rep5.R"

Rmarkdown

---

Split page into three columns (displays best in browser).
R code is in Rmd file.

# r plot code
require(ggplot2)
ggplot(mtcars, aes(x = mpg)) + geom_histogram(fill = "skyblue", alpha = 0.5) + theme_classic()

Praise the lord, I was born to travel
Feeling like Slash in front of the chapel
I'm leaned back with the Les Paul
Shit I smoke is like cholesterol
Spilled dressin' on the vest at the festival
The best of all, had a midget Puerto Rican at my beckon call

Pump the bass in the trunk
It rattled like a baby hand
Except this toy cost 80 grand
And I'm crazy tan, from all the places that I've been
Just from writing words with a pen

---

Inserting text within indented equations
\[ p(x) = \theta^{x} (1 - \theta)^{1-x} ~~\mbox{ for }~~x = 0,1 \]

plotly

HTML widget with plotly and crosstalk

require(pacman)
p_load(plotly, tidyr, crosstalk)

m <- gather(mpg, variable, value, -c(year, cyl))
msd <- highlight_key(m, ~variable)
gg <- ggplot(m, aes(factor(year), value)) + geom_jitter(alpha = 0.3) + labs(x = "Year") + theme_classic()

bscols(widths = c(11, rep(5, 2)), filter_select("id", "Select a variable", msd, ~variable, multiple = F), 
    ggplotly(gg, dynamicTicks = "y") %>% layout(margin = list(l = 30)), plot_ly(msd, x = ~jitter(cyl), 
        y = ~value, alpha = ~cyl, linetype = NULL, mode = "markers", hoverinfo = "text", text = ~paste0("Cyl: ", 
            round(cyl), "\n", variable, ": ", value, "\nYear: ", year)) %>% add_markers(alpha = 0.3) %>% 
        layout(xaxis = list(showgrid = F, title = "Cylinder"), yaxis = list(showgrid = F)))

Select a variable

# example 2 with changing output margins to fill browser


require(htmltools)

# title
plotbg <- tags$html(HTML("<body style=\"background-color: black;\"></body>"))

resource_type <- "detritus"
memi_df <- readr::read_csv("https://raw.githubusercontent.com/darwinanddavis/mybio/master/data/memi_df.csv")
memi_df <- data.frame(memi_df)
memi_df %>% str

'data.frame':   3600 obs. of  7 variables:
 $ X1              : num  1 2 3 4 5 6 7 8 9 10 ...
 $ ME_EVENT        : chr  "skip30_0" "skip30_0" "skip30_0" "skip30_0" ...
 $ Cercs           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ControlImpact   : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ControlDay      : num  1 1 1 1 1 1 1 1 1 1 ...
 $ ControlDay_names: chr  "Skip~30" "Skip~30" "Skip~30" "Skip~30" ...
 $ Time            : num  1 2 3 4 5 6 7 8 9 10 ...

# heatmap -----------------------------------------------------------------
require(viridis)
require(ggthemes)
require(plotly)

me_day_vec <- c("skip30", "skip60", "skip90", "skip120")
me_day_names <- as.factor(c("Skip~30", "Skip~60", "Skip~90", "Skip120"))  # char vec for labelling facets
ttl <- ""
subttl = ""
xlab <- "Time (days)"
ylab <- "Control intensity"
# turn names into function for labeller for facets
me_im_names <- c("No control", "50%", "75%", "90%", "95%", "99%")
dens <- memi_df[, "Cercs"]
yy <- memi_df[, "ControlImpact"]
xx <- memi_df[, "Time"]
facet1 <- memi_df[, "ControlDay"]
p <- ggplot(memi_df, aes(x = xx, y = yy, fill = dens)) + geom_tile(colour = "gray", size = 0.01, width = 2, 
    linetype = 0) + scale_fill_viridis(name = "Density", option = "magma")
p <- p + facet_wrap(~ControlDay_names, nrow = length(me_day_names), ncol = 1, drop = F, labeller = label_parsed)  # use for adding facet labels
# p <-p + facet_wrap(facet1, nrow=3, ncol=1, drop= F)
p <- p + scale_y_continuous(breaks = unique(yy), labels = me_im_names, trans = "reverse")
p <- p + scale_x_continuous(breaks = seq(0, max(xx), 30), expand = c(0, 2))
p <- p + geom_segment(aes(x = 152, xend = 152, y = 4, yend = 4), arrow = arrow(length = unit(0.2, "cm")))
p <- p + theme_calc() + theme(text = element_text(size = 18)) + # labs(title= paste0('Density of ',ttl, ' by ',subttl), y=ylab, x=xlab) +
labs(title = paste0("\n", "\n", ttl), y = ylab, x = xlab, size = 3) + theme(plot.title = element_text(vjust = -7)) + 
    theme(legend.position = "bottom", legend.direction = "horizontal") + theme(legend.text = element_text(size = 12)) + 
    theme(plot.background = element_rect(fill = "black")) + ggpubr::theme_transparent()
# plot_it_gg('black','white')

m <- list(t = 100, r = 1, b = 1, l = 1, padding = 4)

p <- ggplotly(p)

require(htmlwidgets)
h <- p %>% layout(plot_bgcolor = "black", paper_bgcolor = "black", font = list(color = "black"), autosize = T, 
    margin = m) %>% sizingPolicy(padding = 0, browser.fill = TRUE, plotbg)

Crosstalk example 2

# time series plotly
pacman::p_load(dplyr, lubridate, ggplot2, plotly, gridExtra, plyr, ggthemes)
# install.packages('crosstalk')
library(crosstalk)

# load mock data
df <- readr::read_csv("/Users/malishev/Documents/Data/time_series/call_activity/call_activity.csv")
df %>% head

# A tibble: 6 x 4
  Date                Person  Hour Calls
  <dttm>              <chr>  <dbl> <dbl>
1 2018-09-25 00:00:00 Ben        8     1
2 2018-09-26 00:00:00 Rob       16    11
3 2018-09-27 00:00:00 Matt      18    11
4 2018-09-28 00:00:00 Ben       10     8
5 2018-09-29 00:00:00 Rob        9    11
6 2018-09-30 00:00:00 Matt       8     8

xinter <- seq(min(df$Date), max(df$Date), length.out = length(df$Date))

# plot data
p <- ggplot() + geom_vline(mapping = NULL, xintercept = xinter, colour = "grey80", size = 0.03) + geom_point(data = df, 
    aes(Date, Hour, color = Person, size = Calls)) + scale_y_continuous(limits = c(1, 23)) + scale_x_datetime(date_breaks = "1 week", 
    date_minor_breaks = "1 day", date_labels = "%D") + theme(axis.text.x = element_text(angle = 45)) + 
    labs(title = "Calls per hour of day", x = "Date (M/D/Y)", y = "Hour of day") + theme(panel.border = element_blank(), 
    panel.grid.major = element_line(color = "gray"), panel.grid.minor = element_line(color = "light gray"), 
    axis.line = element_line(color = "gray"))
p <- p + theme_hc()
ggplotly(p)

# plotly crosstalk
calls_person <- highlight_key(df, ~Hour)
person_person <- highlight_key(df)

pp <- bscols(widths = 12, p1 <- plot_ly(df, x = ~Date, y = ~Hour, color = ~Person, size = ~Calls, type = "scatter", 
    hoverinfo = "text", text = ~paste0("Date: ", Date, "\nName: ", Person, "\nCalls: ", Calls)) %>% layout(title = "Calls per hour of day", 
    xaxis = list(tickangle = 45, showgrid = T), yaxis = list(range = c(0, 23), showgrid = T), margin = list(l = 0.5)), 
    filter_select("id", "Select hour of day", calls_person, ~Hour, multiple = F), p2 <- plot_ly(calls_person, 
        x = ~Person, color = ~Person, type = "histogram") %>% layout(title = "Calls per person", yaxis = list(showgrid = F)))

pp <- htmltools::tagList(list(p1, p2))
`?`(`?`(tagList))

Sys.setenv(plotly_username = "malishev")
Sys.setenv(plotly_api_key = "uApW9Ar4GpjbEbagDeAn")

ff <- plotly::api_create(p1, username = "malishev")

subplot(p1, p2, nrows = 2)

htmltools::knit_print.shiny.tag.list(pp)

Gather/melt dfs to make dfs plotly friendly

require(tidyr)
require(plotly)
sm <- as.data.frame(EuStockMarkets) %>% gather(index, price) %>% mutate(time = rep(EuStockMarkets %>% 
    time(), 4))

sm %>% plot_ly(x = sm$time, y = sm$price, color = sm$index)

Time

Get just HMS portion of POSIX class

require(dplyr)
Sys.time() %>% format(format = "%H:%M:%S")